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Abstract — We present a corrective subcell averaging technique 
that improves on the accuracy of the volume-averaged finite- 
difference time-domain (FDTD) method in the presence of dis- 
persive material interfaces. The method is based on an alternative 
effective-medium formulation that captures field discontinuities 
at interfaces as electric and magnetic surface currents. In 
calculating the spectra of strongly dispersive Mie scatterers we 
demonstrate that the derived FDTD algorithm is both highly 
efficient and able to approximately restore second order accuracy. 



I. Introduction 

Half a century after its invention by Kane Yee (H the 
finite-difference time-domain (FDTD) method remains a pop- 
ular choice for simulating the propagation of electromagnetic 
waves and their interaction with electronic media JT], J2). 
The simplicity of the algorithm and its low computational 
footprint are contrasted by the use of non-conformal grids, 
which, if field discontinuities are not properly accounted for, 
reduce accuracy from second to first order [0], 02]. This 
not only negates the advantage of the staggered grid Yee- 
algorithm but also severely impacts on the computational cost 
when modeling systems that exhibit geometric features on sub- 
wavelength scales. 

The problem of restoring accuracy of the FDTD scheme 
in the presence of interfaces was first studied in the mi- 
crowave regime 0, J6). Since then a variety of effective- 
permittivity (EP) models have been suggested for the treat- 
ment of field discontinuities at material interfaces, which 
can broadly be classified as either counter-path (CP) or 
volume-polarized (VP) models Q, 0, 0, 0. Of particular 
importance for this work is the VP model proposed by 
Farjadpour et al iflOl . Based on the continuity of the E| 
and Dj^ field components the volume-integrated permittivity 
tensor is derived as e^ 1 = (e^-)P + (£oo) _1 (l — P), where 
P = n <E> n performs a vector-projection onto the face- 
normal of the interface. The application of this non-diagonal 
permittivity tensor requires interpolation of the Yee-centered 
D field to the cell-center and subsequent interpolation of the 
cell-centered E field back onto the Yee-grid, a procedure that 
effectively equates to a smoothing operation with extended 
spatial stencil. Nonetheless, as numerical evidence suggests, 
the spectral accuracy increases to approximately second order, 
reducing the computational cost for problems that involve non- 
dispersive dielectrics (e.g., photonic crystal applications). In 
2007, Deinega et al ifTTl suggested an approach that extends 
this method to the linear dispersive regime. Their algorithm 



uses the decomposition E 



E 



II 



n{E_ 



±.i 



E x ,2) to 



split the electric field into four components, which drive the 
arization currents at the interface. A recent work by Liu 




pola 



Figure 1. (color online) Schematic of the VA+CC FDTD algorithm. E 
and H fields (black) are calculated using the standard VA FDTD algorithm, 
introducing systematic errors at interface-cells due to discontinuities of the 
field. The surface currents 63 j_ and <5K|| (red) correct the errors during 
the electric (left) and magnetic (right) update steps. Calculating the corrective 
currents requires an intermediate step for integrating the surface charge density 
p at the cell center (middle). The algorithm is compatible with the standard 
Yee scheme as the corrections only apply to interface cells. 



et al lfl2ll shows that, in case of dielectric-dispersive media 
boundaries, one does not need to split the normal field as 
long as one can formulate a modified effective polarization 
response for the normal component. The resulting algorithm is 
less general yet computationally more efficient. However, as in 
fTTl . it remains unspecified how the algorithm interfaces with 
the standard Yee-algorithm that could be efficiently employed 
across regions where permittivities are smooth. 

Here, we present an alternative approach (see Fig. Q3 that 
solves the EP curl equations on the Yee-grid using the stan- 
dard volume-averaged (VA) FDTD algorithm but replaces the 
electric field with an approximate field E that is continuous 
across non-dispersive interfaces. The field discontinuities at 
dispersive media interfaces need then to be captured as cor- 
rective electric and magnetic currents S3± and <5K||, which 
are induced by a surface charge field p. Based on this idea 
we first formulate an effective medium theory and then show 
how this EP model translates into a FDTD scheme that offers 
some unique advantages: 1) the algorithm naturally extends 
the standard FDTD scheme by introducing additive current 
corrections that only apply at interface cells; 2) only the 
normal field components are subjected to spatial smoothing 
operations at the interface; and 3) calculating the corrections 
is computationally efficient and requires no alteration of the 
dispersive response functionals (as for example in lfl2l ). In the 
result section, we apply the derived algorithm to the example 
of a highly dispersive Mie scatterer in two-dimensions, demon- 
strating stability and allowing for a comparison of numerical 
errors between the VA+CC (using current corrections), the VA, 
and standard staircasing schemes. 

II. Corrective-Current Subcell Smoothing 

Our starting point are the split-field equations derived by 
Deinega et al [ 11 1 (equations (.T)-(ri) therein! Without loss of 
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generality we write the scalar permittivity as s(lu) — + 
and transform the equations into time-domain. Using a 
slightly different notation, we write 

( £oo )a f E|| - [V x H] N - /iJi[E||] - / 2 J 2 [E||] 
eoo,ifl^i,t=A[VxH]x-J 1 [Bx 4 ] (1) 
Soo,2d t E ±t2 = / a [V x H]x - J 2 [£x )2 ] 

where Jiy 2 = <3tPi/ 2 are functionals of E, describing 
the (isotropic) polarization current response. The symbols 
'_L' and '||' denote vector-projections relative to the interface 
with face-normal n and the notation Jx/ 2 = n • Ji/ 2 is 
introduced for brevity. This formulation of Ampere's law 
implicitly assumes an averaging over a volume-cell that is 
intersected by a boundary between media 1 and 2 with cell- 
filling ratios j\ and / 2 (j\ + / 2 = 1). Angled brackets are 
used throughout this work to denote volume averages of the 

form (goo) = /l£oo,l + /2£oo,2- 

The derivation of ((TJ is straightforward but their translation 
into an efficient and stable finite-difference scheme is not. To 
retain second order accuracy, the field components in the curl 
expression V x H should be calculated on the Yee-grid while 
the projections onto parallel and normal projections require 
interpolation to the cell-center. After calculating the updates 
of the E||, E± t i and E± t 2 components at the cell-center 
the E-field thus needs to be reconstructed and redistributed 
onto the Yee-grid. However, a direct implementation proves 
impractical for the following reason. The cell-centered four- 
field representation and the extended spatial stencil (due to 
interpolation between the grids) fundamentally differs from 
the standard Yee-algorithm. As a consequence the algorithm 
is deployed across the whole grid irrespective of whether cells 
are intersected by media-boundaries or not. This introduces 
unnecessary smoothing operations across the whole grid, in- 
creases the computational cost and requires re-implementing 
the infrastructure typically associated with FDTD frameworks 
(e.g., total-field scattered-field injection, boundary conditions 
etc). 

We therefore seek to derive an alternative formulation where 
the standard Yee scheme can be efficiently applied across 
the domain augmented by corrections that only apply to the 
comparably small number of interface cells. The basis for this 
corrective method is a reformulation of ((T). In introducing new 
variables for the normal electric field and the density of the 
induced surface charges, 

Ei_ = Ej_,i + Ej_, 2 

P = f2^oo,\E_\_,l — /l£oo,2-Ex,2 

Equation (Q~|i can be cast into the form 

( £oo )9 t E|| = [V x H] N - AJi[E,|] - / 2 J 2 [E||] 
(e^ 1 )- 1 dtE ± = [V x H] ± - CiME±,i] - C 2 J 2 [£±, 2 ] (3) 
dtp = hME± t2 ] - f 2 Jx[E ±!l ) 

with Ci/ 2 - (e^r'E' 4/ 2 - The fact that (e^ 1 ± {e^ 1 } 
makes it impossible to reconstruct Ampere's law for the 
discretized (i.e., volume averaged) electric field by adding the 
first two equations. However, we can define an approximate 



electric field 

E = E|, + (e 00 )- 1 <£ M 1 )- 1 n£ ± (4) 

which, in the absence of dispersive currents, is continuous 
across material interfaces and matches E at non-interface cells. 
Combining the first two equations of (0 in this fashion yields 

( £oo )3 f E = VxH-(J[E])-iJ 1 (5) 

We note that apart from the extra current term 5J± we 
now have recovered the volume-averaged curl equation for 
the electric field. The correction i5J^ = n<5Jx compensates 
the error that arises from using volume-averaged permittivities 
and current densities for the normal components. Assuming an 
isotropic response one obtains after some algebra 

SJ± = -/iJi[E] - / 2 J 2 [E] + CiJi[E ±il ] + C 2 J 2 [£x, 2 ] (6) 

for the surface current correction. Its calculation requires the 
scalar fields £?x,i/ 2 that are obtained by projection 

Sx > i = /ie- 1 >1 «s c)n-E + / 1 - 1 C2p) 
£x,2 = / 2 £-y<£oo)n • E - f 2 -\ lP ) 
Inserting these relations into © yields 

5J± =/i ( J? [E, p] — J, [E] ) + h (A [E, p] - J 2 [E] ) (8) 
where we defined 

Jr /2 [E, P ] =Ci/2e M 1 , 1/2 ((eoo)Ji/2[E]±/r/ 2 C 2 /iJi/ 2 [p]) (9) 

This implies that the electric current correction can be cal- 
culated from the currents induced by E and p. The terms in 
^ proportional to Jj/afE] are the volume-averaged normal 
currents, which need to be subtracted from eq. (|5) before 
adding the correct J*/ 2 [E,p] contributions. Applying (0 to 
the equation for the charge field p [see OJ] gives 

In order to complete the update of the magnetic field the 
correct electric field E needs to be recovered from E. This is 
achieved by introducing a corrective magnetic current density 

<5Kii = V x nSE^ 

, (ID 
= Vxn((e^ 1 )(e 00 )-l)n-E 

to Faraday's law 

d t U = -V x E - 5K|| (12) 

rendering the overall scheme bi-isotropic. This completes 
our reformulation of the effective cell-averaged Maxwell's 
equations. The curl equations (O, ( fT2l together with the 
electric and magnetic current corrections (|8} and (Q~Tj and 
the surface charge equation ([Tot form a closed set of equa- 
tions. We achieved our goal of finding an effective medium 
formulation where corrective current densities 63 ± and SKu 
apply at interface cells and vanish whenever permittivities 
vary smoothly across cells. The magnetic current correction 
(5K|| accounts for field discontinuities caused by a jump in 
the static permittivity across the interface, while the electric 
current correction S3± captures all discontinuities induced 
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Figure 2. (color online) Pictorial representation of the action of the C, P 
and Y operators. C interpolates the vector of Yee-centered components F to 
the cell-centered vector F* (left), P projects a cell-centered vector onto the 
face normal (center), and Y interpolates a cell-centered vector back onto the 
Yee-grid (right). 

by the dispersive material response. Notably, calculating the 
induced corrections requires only three additional physical 
fields, namely the interface charge field p and the associated 
induced normal currents J\/2[p\- 

III. Yee-compatible Corrective-Current FDTD 

SCHEME 

We now proceed to translate the equations derived in the 
previous section into a versatile and efficient FDTD scheme. In 
compliancy with the standard Yee-scheme we integrate (|5) and 
( fl2] i in two distinct half-steps by first performing the electric 
field update 

fn+l/2 = £m-l/2 + At(e 0o )- 1 V X H n 

(13) 

-Ai( £oo )- 1 ((J"[E"- 1 / 2 ]}+<5Jl) 
and then the magnetic field update 

jjn+l =H n _ At (y x jj.n+1/2) + Ai<5K™ +1/2 (14) 

To shorten the notation, we implicitly assume that E"+V2 and 
H" are 37V-dimensional vectors (N being the number of Yee- 
cells) aggregating the electric and magnetic field components 
on the staggered subgrids across the problem domain. In this 
formulation the curl-operator V x is a matrix that performs a 
stencil operation at each point of either the electric or magnetic 
subgrid. The permittivity tensor (Eoo) 1 is a diagonal matrix 
that can be precalculated by volume-averaging the permittivi- 
ties at the various positions of the Yee-cube. In a similar way 
(J"[E™ -1 / 2 ]} can be obtained by weighting the contributing 
current vectors J™ with the matrix of precalculated cell- 
filling factors fi. It is important to note that the treatment 
of dispersive currents requires a preceding evaluation of the 
response functionals J"[. . . ], by either integrating appropriate 
auxiliary differential equations (e.g., for the Lorentz pole) |2| 
or by using the piecewise linear recursive convolution (PLRC) 
method lfl3l . 

Following the arguments laid out in the previous section it 
is clear that the corrective currents (5J" and 5K" +1 ^ 2 vanish 
whenever the material constants vary smoothly across cells. 
For these volume cells E — > E and the update equations 
reduce themselves to the dispersive VA FDTD method, which, 
as (eoo} 1 is diagonal can be efficiently integrated using the 
standard Yee-scheme. Within interface cells, on the other hand, 
E differs from the electric field E and a corrective step is 



necessary to accurately account for the discontinuity of the 
normal field component. As shown before the discontinuity in 
the normal component is directly proportional to the surface 
charge density p induced at the interface. Discretizing (fTOb 
results in an update equation for p 

p n+l/2 = p n-l/2 + A^/^-l J*« _ ^ J*™) (1 5 ) 

that requires evaluation of the currents Jj*™ 2 according to ©. 
In difference to the electromagnetic field components, which 
are evaluated on the Yee-grid, p is a cell-centered quantity. We 
therefore need to introduce operators to interpolate between 
the Yee- and cell-centered grids. Figure [2] illustrates the action 
of the Y and C interpolation operators (left and right panel) 
together with the projection operator P (center panel). Applied 
to write (O this yields 

Jl/2 = Cl/ 2 Cl/2( PC ( £ oo)Jv 2 [E"- 1 / 2 ] 

±/r/ 2 C2/iJr/2[p^ 1/2 ]) 

This expression recycles the previously calculated J™^ 2 
currents on the Yee-grid but introduces a charge-current 
J™ /2 [p n ~ 1/2 } that, using the same current-functional, is eval- 
uated at the cell center. To improve smoothness of the fields 
under the projection/interpolation operation we multiply J™, 2 
with the (eoo) tensor, which is already available on the Yee- 
grid. In contrast, the coefficients £1/2, £ ^\/2 an< ^ fi/2 an< ^ 
the face-normal n are parameters that are defined at the cell- 
center (see Fig. [8j. As (TTol i can be evaluated on-the-fly, the 
only additional physical fields that need to be stored at the 
cell-center are p and its induced currents J?, 2 [p n_1 ' 2 ]. 

With the surface charge and its currents known, it becomes 
possible to compute the corrections <5J™ and <5Kjj +1 ^ 2 that 
enter the update equations ( fT3] > and (TT4l . However, the order 
of operators (and hence the discretization) is ambiguous, and, 
as the scheme is corrective, can impact on the stability of the 
scheme. A numerical analysis of the computational errors sug- 
gests that J™^ 2 is best multiplied with the (£oo) tensor before 
centering to the grid. This is due to the fact that the normal 
component of (eoo)E retains smoothness across adjacent cells 
with different e^. Further, to maintain consistency between 
the Yee and cell-centered update equations (I131 l. (IT4b and ([13] ) 
we assign parameters as indicated by Fig. [8] This allows us 
to write 

(e^SJl = - (e 00 )- 2 (/iYPC{e 00 )J' 1 i - ^YPC^))^ 1 
+ Ynf 1 J* n +Yn/ 2 J* n 

(17) 

where volume filling factors in the first line are applied after 
centering onto the Yee-grid, and, for the second line, directly 
at the cell-center. 

The discretization of the magnetic current requires both 
terms in (fTTT l to be interpolated to the center before spreading 
them out again onto the Yee-grid. We obtain 

<5Kf] +1/2 = Vx^ +1/2 (18) 

with 

5& +x/2 = (Y(£ _ 1} (eoo) -i Y )PC( eoo )E" +1 / 2 (19) 
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Algorithm 1 Sequence of field updates (VA+CC) 

• n+1/2 (on Yee-grid): 

- VA: update E™" 1 / 2 _> E n+1 / 2 w /o [(E)] 

- CC: add correction [Q3] 

• n+1/2 (on centered- grid): 

- CC: update p™" 1 / 2 ->• p n + 1 /2 [gjj] 

- CC: evaluate J™ +1 [/>" +1 / 2 ] (for next cycle) 

• n+1 (on Yee-grid) 

- VA: evaluate J" +1 [E' i+1 / 2 ] (for next cycle) 

" ten 



(a) 



- VA: update H" ->• H n+1 w/o <5K 

- CC: add correction <SK™ +1/2 [([TD] 
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Fundamentally, both the electric and magnetic current cor- 
rections can be calculated on-the-fly. As the corrections only 
apply to interface cells, they can be added in a separate step 
to the update equations. This means that the update equations 
of the VA FDTD scheme can be deployed across the whole 
grid, followed by oversampling steps that perform the current- 
corrections (CC) for interface cells only. The complete update 
sequence for the VA+CC algorithm is shown in Alg. Q] 

As each step can be associated with a loop over cells, 
it becomes evident that the current-correction (CC) steps 
augment those related to the VA FDTD scheme. As the CC 
steps only apply to interface cells, the computational overhead 
of the VA+CC FDTD scheme is not significant unless the 
number of interface cells becomes comparable to the number 
of volume cells. 

IV. Results 

To verify the accuracy of our method we compare our 
numerical calculations with the Mie scattering cross section 
of an infinitely extended strongly dispersive cylinder excited 
by a TM plane-wave. The dielectric function describing the 
response of the cylinder consists of a single Lorentzian reso- 



nance at Ar = 0.25i? and a background dielectric constant 
of £oo = 4, where R is the radius of the cylinder. Figure [3^ 
shows real and imaginary parts of the complex permittivity 



e(\- 1 ) = e'(\- 1 )+is"(\- 1 ) 



+ 2.5Aq (Ao - A 



i0.057r _1 A -1 ) -1 together with the analytically calculated 
scattering cross-sections for scatterers with and without the 
dispersive contribution x(A _1 ) (Fig. [3J)). 

The numerical setup of the 2D calculation is depicted in Fig. 
|U An incident field is injected into the system at one of the 
boundaries of a standard Total-Field-Scattered-Field (TFSF) 
box [2|. The energy flux E x H of the scattered field is then 
recorded at the boundary of a box located outside of the TFSF 
box. The computational region is terminated with perfectly 
matched layers (PML) [2] which nearly completely attenuate 
any reflections caused by the computational boundary. After 
the simulation, the scattering spectrum can be recovered by 
Fourier-transforming the recorded energy flux. 

Figure (top) shows the difference between the analytic 
and numerical scattering cross sections obtained by numerical 
simulation with a resolution of eight Yee-cells per cylinder 
radius. The results of the VA+CC FDTD scheme (dashed red 
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Figure 3. (color online) (a) the real (dashed) and imaginary (dotted) part of 
the dielectric function of an infinitely long cylinder and (b) the analytically 
calculated scattering cross section (solid). The dotted curve in (b) represents 
the scattering cross section of a cylinder with a purely static dielectric constant 
of £00 = 4 (indicated by the thin dotted line in (a) 




Figure 4. Computational setup: an incident pulse is injected on the inner left 
boundary of a TFSF box. The pulse interacts with the scatterer and leaves the 
TFSF box outside of which the energy flux of the scattered field is recorded 
(dotted line marked "DIAG"). The boundary of the computational region is 
terminated with perfectly matched layers (PML). The current corrections (CC) 
<5J^ and 5Kn are only applied at the surface of the cylinder (thick dashed 
line). 



line) are in better agreement with the analytical calculation 
than the VA FDTD scheme (dotted green line) throughout the 
spectrum. For comparison, the result of a simple staircased 
FDTD scheme was included in the figure (dash-dotted blue 
line). By selectively disabling either the current correction 
(5J^ or (5K|| and subtracting the result from the VA FDTD 
scheme, the contributions of the charge corrections to the 
spectrum were quantified (Fig. ^>). The contribution of <5Kn 
is significant in the entire recorded spectral range, while, the 
contribution of <5Jj_ shows a peak at a frequency which is 
slightly offset to the resonance frequency of the Lorentzian 
(indicated by the vertical dotted line). To illustrate the spatial 
dependence of the corrections and the charge density we plot 
contour images of the charge field p (Fig. [(J), the energy 
density of the electric correction 83 ± ■ E (Fig. [(J), and the 
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Figure 5. (color online) (a) difference between the analytic and numerical 
scattering cross section spectrum of an infinitely long dispersive Mie cylinder 
calculated with a VA FDTD scheme with (red dashed) and without (green 
dotted) charge corrections. The result of a staircased FDTD scheme is shown 
for reference (dashed-dotted blue). The thin dotted black horizontal line is a 
guide to the eye. (b) the contribution of the individual charge corrections to 
the cross section spectrum 
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Figure 6. (color online) snapshot of (a) the volume filling factor /2 of the 
interface cells (b) the charge field p (c) the electric current correction energy 
<5J x ■ E (d) the magnetic current correction energy <JKii ■ H. 

energy density of the magnetic correction <5Kii ■ H (Fig. |6}i). 
Whereas the corrections associated with the charge density 
and electric field correction are stored at the cell center, the 
correction associated with 6~K\\ ■ H is calculated from Yee- 
centered quantities and therefore appears to be smeared out 
over several adjacent cells. 

To investigate the convergence behavior of the charge 
correction algorithm, numerical simulations with increasing 
resolution N were conducted. The RMS error for each sim- 
ulation was obtained, by comparing the numerical scattering 
cross section spectrum with the analytical result (Fig.|7^). The 
VA (green diamonds) and the staircasing (blue circles) FDTD 
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Figure 7. (color online) RMS error vs. increasing pixels per cylinder radius 
TV. (a) the RMS error of staircasing (blue circles), VA (green diamonds) 
and VA+CC FDTD (red squares) schemes, (b) the relative change in error 
when enabling either the 83 j_ (cyan circles) or <5Ki i (black diamonds) current 
correction compared to the VA FDTD scheme error. 



scheme produce errors that are significantly larger than those 
of the VA+CC scheme (red squares), which shows a decrease 
of error proportional to oc N~ 2,1 , This means that for this 
particular system VA+CC achieves an accuracy that exceeds 
second order. 

The overall error reduction is achieved by the combined ac- 
tion of the corrections 53 ± and SK.n. To quantify their impact, 
we isolate their contributions to the error reduction as shown 
in the lower graph of Fig. [7] For low resolutions (smaller than 
s» 10 pixels per cylinder radius) the charge correction <5Km 
is dominant in reducing the error. In this resolution regime, 
volume averaging of the static dielectric constant (£oo) for the 
perpendicular electric field component Ej_ is the source of 
the large error in the VA scheme. The charge correction <5Kh 
corrects this error. With higher resolutions, however, the error 
of the volume averaged dynamic response (J) becomes larger 
and needs to be compensated by 5J±. 

Finally, we compare the computational cost (memory and 
processing time) for the different schemes. The results are 
summarized in Fig. [8] The staircasing scheme only requires 
the static epsilon at each Yee-cell position of the E-field 
and the three vectorial fields E, H, Ji. The VA algorithm 
additionally stores the filling factors /i at each Yee-cell 
position of the E-field. The VA+CC scheme is identical to 
the VA scheme for non-interface cells requiring 15 scalar 
components. At interface cells the VA+CC scheme requires 
an additional 7 scalar components for storing p, J\ [p] , n, /i 
and £00,1/2- The comparision of computation time indicates 
an almost identical performance for the staircase and VA 
schemes. VA+CC delivers the same performance for volume 
cells but requires additional computational steps for interface 
cells, resulting in w 50% overhead in the per cell processing 
time. A 50% increase in per cell memory and processing time 
seems significant but rarely matters for practical applications 
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Figure 8. Computational cost (memory and CPU) of staircasing, VA and 
VA+CC FDTD algorithms. The various fields and parameters are either 
assigned to cell-centered (CG) positions or to the Yee-grid (YG). VA+CC 
requires same storage as VA for volume cells but carries an overhead of 
-50% for interface cells. The per-cell storage values are given in QWORDs, 
the per-cell time in microseconds. 



as the surface to volume ratio is typically small. For the Mie 
scattering simulations presented in Fig.[5]fbr example (8 cells 
per radius) the increase in computation time of the VA+CC 
algorithm is < 1% (compared to VA) as the interface/volume 
cell ratio is m 0.6%. 
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V. Conclusion 

In summary we presented an effective-medium theory that 
takes a corrective approach to the cell-averaged Maxwell's curl 
equations. The theory holds for static and linear dispersive 
permittivities and captures the field discontinuities inside a 
cell in form of surface current corrections, which can be 
calculated by integrating a surface charge equation alongside 
the volume-averaged curl equations. We derived a computa- 
tionally efficient FDTD algorithm that allows deploying a Yee- 
centered stencil across the domain followed by surface current 
corrections that apply to interface cells. The improvement in 
accuracy is quantified by calculating spectral scattering cross- 
sections of strongly dispersive Mie scatterers. The extracted 
error exponents indicate that the algorithm approximately 
restores second order accuracy. The work presented is relevant 
in the current context of nano-photonic research and may pave 
the way to the development of novel pertubative techniques for 
solving Maxwell's equations. 
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